---
title: "Online appendix"
subtitle: "Breaking out of Legacy Mobilization Networks: How the Internet Reaches and
Activates the Politically Disengaged"
author: "Francesco Bailo (francesco.bailo@sydney.edu.au)"
institution: "The University of Sydney"
geometry: margin=1.5cm
output:
  bookdown::pdf_document2:   
    toc: false
    keep_tex:  true
---

```{r setup, include=FALSE}

knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE, cache = FALSE)

knitr_cache_path = 'main_cache/latex/'

library(tidyverse)
library(sf)
library(knitr)
library(kableExtra)
library(ggrepel)
library(ggpp)
library(mcp)
library(cowplot)
library(scales)
library(modelsummary)
library(CARBayes)
library(coda)

ita_cities.df <- 
  read.csv("data/ita-cities/it.csv") %>%
  dplyr::filter(city %in% c("Milan", "Turin", "Bologna", "Padova", "Genoa", 
                            "Florence",
                            "Rome", "Naples", "Bari", 
                            "Palermo", "Catania", "Messina",
                            "Cagliari"))

```

```{=tex}

\renewcommand{\thetable}{A\arabic{table}}
\renewcommand{\thefigure}{A\arabic{figure}}

```

**Cite**: Bailo, F. (2025). Breaking out of legacy mobilization networks: How the Internet reaches and activates the politically disengaged. **Journal of Information Technology & Politics**. [https://doi.org/10.1080/19331681.2025.2494116](https://doi.org/10.1080/19331681.2025.2494116)

Replication materials are available on the Harvard Dataverse [https://doi.org/10.7910/DVN/MNHJTQ](https://doi.org/10.7910/DVN/MNHJTQ).


# Definition of the hexagonal grid

A hexagonal cell size of approximately 15 km is defined with `italy.sf` as the bounding box. 

```{r echo = T, eval = F}

italy_hex.sf <-
  sf::st_make_grid(italy.sf,
                   15000,
                   what = "polygons",
                   square = FALSE)

italy_hex.sf <- 
  sf::st_sf(index = 1:length(lengths(italy_hex.sf)), 
            italy_hex.sf)


```

The resulting `italy_hex.sf` is then cropped using `italy.sf` as mask. 

```{r echo = T, eval = F}

italy_hex_cropped.sf <-
  sf::st_intersection(italy_hex.sf, 
                      italy.sf %>% st_make_valid())

```

```{r}

load("data/out/define-hex.RData")

```


```{r hex-grid, fig.cap=sprintf("Hexagonal cells of size 15km (n=%s)", nrow(italy_hex_cropped.sf)), fig.height=10, out.height="35%", fig.align = "center", eval = F}

The resulting grid is shown in Figure \@ref(fig:hex-grid).

ggplot(italy_hex_cropped.sf) +
  geom_sf()

```

# Intersection between census sections and hexagonal cells for count variables

```{r}

capture.output(
  knitr::load_cache(label = "sez2011-lonlat.sf",
                    object = 'sez2011_lonlat.sf',
                    path = knitr_cache_path),
  file='NUL')

```


To map census count variables measured at the level of census sections with hexagonal cells, I proceeded as follows:

1. I computed the centroid, or its geographic centre, of the census sections;
2. I intersected centroids and hexagonal cells;
3. For each hexagonal cell, I summed the count variables pertaining to the census sections with a centroid contained within the cell.

![](figure/census_to_hex_grid.pdf)

* Sections of the 2011 census: `r format(nrow(sez2011_lonlat.sf), big.mark=",")`;
* Hexagon cells associated with at least one census section: `r length(unique(sez2011_lonlat.sf$hex))`;
* Distribution of census sections by hexagon cells: Min `r paste0(names(summary(as.numeric(table(sez2011_lonlat.sf$hex)))), ": ", summary(as.numeric(table(sez2011_lonlat.sf$hex))))`


```{r}

regions.sf <- 
  sf::read_sf("data/istat/Reg01012020_g/Reg01012020_g_WGS84.shp",
              quiet = TRUE)

capture.output(
  knitr::load_cache(label = "hex-ranking-on-pop-density",
                    object = 'italy_hex_cropped.sf',
                    path = knitr_cache_path),
  file='NUL')

```

# Intersection between municipality and hexagonal cells for vote variables

To map municipality *comune* count variables onto hexagon cells, I proceeded as follows: 

1. For each municipality, vote counts were distributed across all the census sections found within the municipality boundaries proportionally to the size of the population residing in each cell.

2. Vote counts estimated for each census section were then assigned to hexagonal cells containing the sections' centroid, with the same method used for the census variables.


# Cell selection for the spatial analysis

Population density varies significantly across the country and cells. I created an ordinal variable to capture this variation with a 6-point scale index. The variables with key statistics for each group are presented in Table \@ref(tab:pop-density-hex).  

```{r pop-density-hex}

italy_hex_cropped.sf$P1[is.na(italy_hex_cropped.sf$P1)] <- 0 

italy_hex_cropped.df <- 
  italy_hex_cropped.sf

sf::st_geometry(italy_hex_cropped.df) <- NULL

italy_hex_cropped.df$density <- 
  italy_hex_cropped.df$P1 / (as.numeric(max(italy_hex_cropped.df$area,na.rm=T)) * 1e-6)

italy_hex_cropped.df %>%
  dplyr::mutate(`ordinal` = str_extract(P1_cut_6, "^[1-6]")) %>%
  dplyr::group_by(`ordinal`) %>%
  dplyr::summarise(cells = n(),
                   `tot` = sum(P1),
                   `%` = round(sum(P1) / sum(italy_hex_cropped.df$P1) * 100, 2),
                   `min pop` = min(P1),
                   `max pop` = max(P1),
                   `min den` = round(min(density),2),
                   `max den` = round(max(density),2),
                   `mean den` = round(mean(density),2),
                   `median den` = round(median(density),2),) %>%
  kbl(booktabs = T, linesep = "", 
      col.names = c("ordinal", "n. cells",
                    "tot", "%",
                    "min", "max",
                    "min", "max","mean","median"),
      caption = "Population and density of hexagon cells") %>%
  add_header_above(c(" " = 2, "population" = 4, " density (per km2)" = 4)) %>% 
  kable_styling(latex_options = "striped", font_size = 7)

```

The 864 cells (46%) with an ordinal value of 1 on the density index are excluded from the spatial analysis. These  cells cumulatively contain about 7.5% of the national population (or 4,454,720 people) but are considered sparsely populated as their median density is only 24.4 people per km$^2$. 

The geographic distribution based on the population density ordinal variable is mapped in Figure \@ref(fig:plot-hex-pop-cuts).

```{r plot-hex-pop-cuts, fig.cap = "Population density ordinal variable", fig.height=10, out.height="35%", fig.align = "center"}

mycol <- 
  c("#ca0020", "#ffffcc", "#a1dab4", "#41b6c4", "#2c7fb8", "#253494")

ggplot() +
  geom_sf(data =
            italy_hex_cropped.sf %>%
            sf::st_transform(4326) %>%
            dplyr::mutate(ordinal = str_extract(P1_cut_6, "^[1-6]")),
          aes(fill = ordinal), colour = NA) +
  scale_fill_manual(values = mycol) +
  geom_sf(data = regions.sf %>%
            sf::st_transform(4326), 
          fill = NA, size = .5) +
  geom_point(data = ita_cities.df, aes(lng, lat), shape = 1, size = 6,
             alpha = .35) +
  geom_text_repel(data = ita_cities.df, aes(lng, lat, label = city), 
                  fontface = "bold", 
                  nudge_x = c(-.5, -.5, -.5, -.5, -.5, 0,
                              2, -2, .5, 1.5, -.5, 1.5, .2), 
                  nudge_y = c(-.5, 1.5, -.5, 1.5, .5, -.5,
                              .2, -.5, .5, -.05, .8, -.2, -.5),
                  max.overlaps = 15, alpha = .75,
                  point.padding = .5) +
  labs(fill = NULL) +
  theme(legend.position = 'bottom',
        legend.text=element_text(size=8))


```


```{r eval = F}

cowplot::plot_grid(
  vars_1.df %>%
    dplyr::mutate(tot_pop = sum(P1)) %>%
    dplyr::group_by(P1_cut_6) %>%
    dplyr::mutate(P1_cut_6_desc = sprintf("%s,  pop: %s%%", P1_cut_6, 
                                          round(sum(P1) / tot_pop[1] * 100))) %>%
    dplyr::filter(P1_cut_6 != "6") %>%
    ggplot(aes(x = sqrt(NP_949_perc), y = sqrt(rsvp_sum_2007_2012_perc))) + geom_point() + geom_smooth(method = 'lm') + 
    facet_wrap('P1_cut_6_desc', scales = "free") +
    theme_bw() +
    stat_cor(method = "pearson") +
    labs(x = "staff not-for-profit ass., %", y = "RSVPs meetup, %"),
  vars_1.df %>%
    dplyr::mutate(tot_pop = sum(P1)) %>%
    dplyr::group_by(P1_cut_6) %>%
    dplyr::mutate(P1_cut_6_desc = sprintf("%s,  pop: %s%%", P1_cut_6, 
                                          round(sum(P1) / tot_pop[1] * 100))) %>%
    dplyr::filter(P1_cut_6 != "6") %>%
    ggplot(aes(x = sqrt(NP_949_perc), y = sqrt(rsvp_sum_2012_2013_perc))) + geom_point() + geom_smooth(method = 'lm') + 
    facet_wrap('P1_cut_6_desc', scales = "free") +
    theme_bw() +
    stat_cor(method = "pearson") +
    labs(x = "staff not-for-profit ass., %", y = "RSVPs meetup, %"),
  labels = c("2007-12", "2012-13"),
  label_size = 12,
  label_x = 0, label_y = 0,
  hjust = -0.5, vjust = -0.5)
```

# Meetup event data

Meetup event data were directly collected from the API of Meetup.com at different times between 2014 and 2018 using search terms related to the Movement (`beppe+grillo`, `movimento+5+stelle`, `movimento+cinque+stelle`, and `m5s`)  

```{r load-meetup}

load("data/meetup/meetup_event.RData") 

load("data/meetup/meetup_venue.RData") 

load("data/meetup/meetup_group.RData") 

load("data/meetup/meetup_venue_googleapi.RData")

load("data/out/load-intersect-meetup.RData")

meetup_events.df$lon <-
  st_coordinates(meetup_venue.sf)[,1][match(meetup_events.df$venue,
                                            meetup_venue.sf$venue_id)]

meetup_events.df$lat <-
  st_coordinates(meetup_venue.sf)[,2][match(meetup_events.df$venue,
                                            meetup_venue.sf$venue_id)]

```

```{r meetup-geo-event, fig.cap = "Geocoded meetup events present in the data", fig.height = 10, dev = 'jpeg', dpi = 300, out.height="35%", fig.align = "center"}

ggplot() +
  geom_point(data = meetup_events.df %>%
               dplyr::group_by(group_id) %>%
               dplyr::arrange(date) %>%
               tidyr::fill(lon, lat, .direction = "downup"),
             aes(x = lon, y = lat, size = yes_rsvp_count),
             alpha = .55, colour = '#1a9850') + 
  scale_size_area(max_size=3) + 
  geom_sf(data = regions.sf %>%
            sf::st_transform(4326), 
          fill = NA, size = .5) +
  geom_point(data = ita_cities.df, aes(lng, lat), shape = 1, size = 6,
             alpha = .85) +
  geom_text_repel(data = ita_cities.df, aes(lng, lat, label = city), 
                  fontface = "bold", 
                  nudge_x = c(-.5, -.5, -.5, -.5, -.5, 0,
                              2, -2, .5, 1.5, -.5, 1.5, .2), 
                  nudge_y = c(-.5, 1.5, -.5, 1.5, .5, -.5,
                              .2, -.5, .5, -.05, .8, -.2, -.5),
                  max.overlaps = 15, alpha = .75,
                  point.padding = .5) +
  coord_sf(xlim = c(6.7499552751, 18.4802470232), 
           ylim = c(36.619987291, 47.1153931748)) +
  labs(size = "event's RSVPs") +
  theme(legend.position = 'bottom')

```

* Number of meetup unique events: `r format(nrow(meetup_event), big.mark = ",")`;
* Meetup events by data collection date: `r paste0(names(table(as.Date(meetup_event$timestamp))),": ", format(table(as.Date(meetup_event$timestamp)), big.mark = ","))`
* First meetup event (by creation time): `r min(as.POSIXct(meetup_event$created, tz = "UTC"))`;
* Last meetup event (by creation time): `r max(as.POSIXct(meetup_event$created, tz = "UTC"))`;
* Number of venues: `r format(nrow(meetup_venue), big.mark = ",")`;
* Number of meetup groups: `r length(unique(meetup_group$group_id))`;
* Number of meetup groups with at least one event: `r length(unique(meetup_event$group_id))`;
* Number of meetup groups that at least once were found not to be private: `r length(unique(meetup_group$group_id[meetup_group$visibility == "public"]))`.

```{r fig.cap = "Meetup RSVP count by week", fig.width = 7, fig.height = 2}

capture.output(
  knitr::load_cache(label = "change-point-original",
                    object = 'meetup_rsvp_count_week.df',
                    path = knitr_cache_path),
  file='NUL')

ggplot(meetup_rsvp_count_week.df, aes(x = as.Date(week), y = sum_yes_rsvp_count)) +
  geom_bar(stat = 'identity') +
  labs(x = NULL, y = "RSVPs")

```

# Intersection between meetup event data and hexagonal cells for RSVP counts

To associate meetup activity with hexagon cells, I intersect the locations of venues with the hexagonal grid. Events that are organised at a venue located within the boundary of a cell and their RSVPs are associated with that cell (see Figure \@ref(fig:meetup-venues-hex)). 

```{r meetup-venues-hex, fig.cap = "Intersection between meetup event data and hexagonal cells for RSVP counts"}

load("data/out/load-intersect-meetup.RData")

this_bb.v <- 
  italy_hex.sf %>%
  dplyr::filter(index == 595) %>%
  sf::st_buffer(25000) %>%
  sf::st_bbox()

this_hex.sf <-
  st_crop(italy_hex.sf, this_bb.v)

this_venue.sf <-
  st_crop(meetup_venue.sf %>% 
            st_transform(32632), this_bb.v)

res <- 
  st_within(this_venue.sf,
            this_hex.sf)

this_venue.sf$hex <- 
  factor(unlist(lapply(res, function(x) 
    return(ifelse(length(x)==0, NA, x)))))

cowplot::plot_grid(
  ggplot() +
    geom_sf(data = this_venue.sf) +
    geom_sf(data = this_hex.sf %>% 
              dplyr::filter(index != 558), alpha = .75) +
    guides(colour = "none") +
    ggplot2::theme_bw() + 
    ggplot2::theme(axis.text = ggplot2::element_blank(), 
                   axis.ticks = ggplot2::element_blank(), 
                   axis.title = ggplot2::element_blank(), 
                   legend.key = ggplot2::element_blank(), 
                   panel.background = ggplot2::element_rect(fill = "white", 
                                                            colour = NA), 
                   panel.border =  ggplot2::element_blank(), 
                   panel.grid = ggplot2::element_blank()),
  ggplot() +
    geom_sf(data = this_venue.sf) +
    geom_sf(data = this_hex.sf %>% 
              dplyr::filter(index != 595), alpha = .75) +
    guides(colour = "none") +
    ggplot2::theme_bw() + 
    ggplot2::theme(axis.text = ggplot2::element_blank(), 
                   axis.ticks = ggplot2::element_blank(), 
                   axis.title = ggplot2::element_blank(), 
                   legend.key = ggplot2::element_blank(), 
                   panel.background = ggplot2::element_rect(fill = "white", 
                                                            colour = NA), 
                   panel.border =  ggplot2::element_blank(), 
                   panel.grid = ggplot2::element_blank()),
  ggplot() +
    geom_sf(data = this_venue.sf) +
    geom_sf(data = this_hex.sf %>% 
              dplyr::filter(index != 650), alpha = .75) +
    guides(colour = "none") +
    ggplot2::theme_bw() + 
    ggplot2::theme(axis.text = ggplot2::element_blank(), 
                   axis.ticks = ggplot2::element_blank(), 
                   axis.title = ggplot2::element_blank(), 
                   legend.key = ggplot2::element_blank(), 
                   panel.background = ggplot2::element_rect(fill = "white", 
                                                            colour = NA), 
                   panel.border =  ggplot2::element_blank(), 
                   panel.grid = ggplot2::element_blank()),
  ncol = 3)


```

# Intersection between meetup event data and municipality boundaries for RSVP counts

The association between meetup event data and municipalities is obtained by intersecting the boundaries of the municipalities with the hexagonal grid. All the meetup events and corresponding RSVPs organised within the boundaries of the hexagonal cells are then passed over to the associated municipalities. Two geographies intersect if they have any number of points in common. As a result, a municipality can also be associated with events that are not within its boundaries in its close proximity (up to 15 km from the boundaries). Three examples are illustrated in Figure \@ref(fig:meetup-hex-municipality). 

```{r meetup-hex-municipality, fig.cap = "Intersection between meetup event data and municipality boundaries for RSVP counts"}

comuni8092.sf <-
  read_sf("data/istat/Com2011_g/Com2011_g_WGS84.shp")

this_comuni.sf <-
  st_crop(comuni8092.sf, this_bb.v)

res1 <-
  sf::st_intersects( this_comuni.sf %>%
                       dplyr::filter(COMUNE == "Reggio nell'Emilia"),
                     this_hex.sf)

res2 <-
  sf::st_intersects( this_comuni.sf %>%
                       dplyr::filter(COMUNE == "Castellarano"),
                     this_hex.sf)

res3 <-
  sf::st_intersects( this_comuni.sf %>%
                       dplyr::filter(COMUNE == "Rubiera"),
                     this_hex.sf)

cowplot::plot_grid(
  
  ggplot() +
    geom_sf(data = this_comuni.sf, colour = "gray") +
    geom_sf(data = this_hex.sf, fill = NA) +
    ggplot2::theme_bw() + 
    ggplot2::theme(axis.text = ggplot2::element_blank(), 
                   axis.ticks = ggplot2::element_blank(), 
                   axis.title = ggplot2::element_blank(), 
                   legend.key = ggplot2::element_blank(), 
                   panel.background = ggplot2::element_rect(fill = "white", 
                                                            colour = NA), 
                   panel.border =  ggplot2::element_blank(), 
                   panel.grid = ggplot2::element_blank()),
  
  ggplot() +
    geom_sf(data = this_comuni.sf %>%
              dplyr::filter(COMUNE == "Reggio nell'Emilia"), 
            colour = "black", fill = "white") +
    geom_sf(data = this_hex.sf, fill = NA) +
    geom_sf(data = this_hex.sf %>%
              dplyr::filter(index %in% this_hex.sf$index[res1[[1]]]),
            alpha = .75) +
    ggplot2::theme_bw() + 
    ggplot2::theme(axis.text = ggplot2::element_blank(), 
                   axis.ticks = ggplot2::element_blank(), 
                   axis.title = ggplot2::element_blank(), 
                   legend.key = ggplot2::element_blank(), 
                   panel.background = ggplot2::element_rect(fill = "white", 
                                                            colour = NA), 
                   panel.border =  ggplot2::element_blank(), 
                   panel.grid = ggplot2::element_blank()),
  
  ggplot() +
    geom_sf(data = this_comuni.sf %>%
              dplyr::filter(COMUNE == "Castellarano"), 
            colour = "black", fill = "white") +
    geom_sf(data = this_hex.sf, fill = NA) +
    geom_sf(data = this_hex.sf %>%
              dplyr::filter(index %in% this_hex.sf$index[res2[[1]]]),
            alpha = .75) +
    ggplot2::theme_bw() + 
    ggplot2::theme(axis.text = ggplot2::element_blank(), 
                   axis.ticks = ggplot2::element_blank(), 
                   axis.title = ggplot2::element_blank(), 
                   legend.key = ggplot2::element_blank(), 
                   panel.background = ggplot2::element_rect(fill = "white", 
                                                            colour = NA), 
                   panel.border =  ggplot2::element_blank(), 
                   panel.grid = ggplot2::element_blank()),
  
  ggplot() +
    geom_sf(data = this_comuni.sf %>%
              dplyr::filter(COMUNE == "Rubiera"), 
            colour = "black", fill = "white") +
    geom_sf(data = this_hex.sf, fill = NA) +
    geom_sf(data = this_hex.sf %>%
              dplyr::filter(index %in% this_hex.sf$index[res3[[1]]]),
            alpha = .75) +
    ggplot2::theme_bw() + 
    ggplot2::theme(axis.text = ggplot2::element_blank(), 
                   axis.ticks = ggplot2::element_blank(), 
                   axis.title = ggplot2::element_blank(), 
                   legend.key = ggplot2::element_blank(), 
                   panel.background = ggplot2::element_rect(fill = "white", 
                                                            colour = NA), 
                   panel.border =  ggplot2::element_blank(), 
                   panel.grid = ggplot2::element_blank()),
  
  ncol = 4
  
)



```

# Change point regression model

```{r}

load("data/out/mcp-fit-original.RData")

```

Two change points in the time series of 386 weekly counts of RSVPs to meetup events (`sum_yes_rsvp_count`) organised everywhere between `r format(as.Date(fit$data$week[1]), '%d %B %Y')` and `r format(as.Date(fit$data$week[length(fit$data$week)]), '%d %B %Y')` are identified using the [mcp](https://lindeloev.github.io/mcp/) R package for multiple change point Bayesian regressions and the following model specification:

```{r echo = T,  eval = F}

model <- list(
  sum_yes_rsvp_count ~ 1 + week_num,
  ~ 1 + week_num,
  ~ 1 + week_num 
)

```

Posterior distribution as well as sampling diagnostics are presented below. `Rhat` is the [Gelman-Rubin convergence diagnostic](https://www.rdocumentation.org/packages/coda/versions/0.19-3/topics/gelman.diag) (lack of convergence is indicated by values above 1) while `n.eff` is the [effective sample size](https://mc-stan.org/docs/2_18/reference-manual/effective-sample-size-section.html). 

Figure \@ref(fig:mcm-plot-fit) shows posteriors and convergence of the two change point parameters.

```{r mcm-plot-fit, fig.width = 7, fig.height = 2.5, fig.cap = "Posteriors and convergence of multiple change point analysis"}

plot_pars(fit, regex_pars = "cp_")

```

The top panel of Figure \@ref(fig:change-points) include the count of weekly RSVPs as points, fitted lines and the two change points posterior density for each of the three chains. The bottom panel of Figure \@ref(fig:change-points) presents the same data points of weekly counts as histograms and the two change points of `r format(week_seq[round(cp_1)], '%d %B %Y')` and `r format(week_seq[round(cp_2)], '%d %B %Y')`.

```{r change-points, fig.cap = "Estimation of change points in the time series of RSVPs weekly counts", fig.width = 10, fig.height = 5}

capture.output(
  knitr::load_cache(label = "change-point-original",
                    object = 'cp_1',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "change-point-original",
                    object = 'cp_2',
                    path = knitr_cache_path),
  file='NUL')

cowplot::plot_grid(
  plot(fit) + 
    theme(axis.title.x = element_blank(),
          axis.text.x = element_blank(),
          axis.ticks.x = element_blank()) +
    labs(y = "RSVPs per week") +
    scale_x_continuous(breaks = c(23, 75, 127, 180, 232, 284, 336)),
  meetup_rsvp_count_week.df %>%
    ggplot(aes(x = as.Date(week), y = sum_yes_rsvp_count)) +
    geom_histogram(stat ='identity') +
    geom_vline(xintercept = week_seq[c(round(cp_1), round(cp_2))], linetype = 2) +
    labs(y = "RSVPs per week", x = NULL) +
    scale_x_date(limits = 
                   as.Date(as.character(
                     c(fit$data$week[1],
                       fit$data$week[length(fit$data$week)])))),
  ncol = 1)

```

# Spatial analysis: regression models

I specified three models to conduct the spatial analysis. 

1. A model predicting meetup RSVPs using two social capital variables: turnout and presence of associations;

2. A model predicting meetup RSVPs using four demographics (as a percentage of the population): population with a university degree, population unemployed, population working as homemakers, and population over 65.

3. A model predicting M5S vote using the same four demographics and meetup RSVPs. In this model all independent variables are scaled to facilitate the visual interpretation of the coefficients. 

An ordinal variable controlling for population density is added to every specification. Additionally, all RSVP variables are subjected to a square transformation assuming a declining marginal effect of additional RSVPs. The variable capturing the difference of RSVPs between the two periods is the difference of two square values.  

All models are fit using the function `S.CARleroux` from the R package CARBayes and the following arguments:

```{r eval = F, echo = T}

S.CARleroux(
  formula,
  data,
  family = "gaussian",
  W = hex.m,
  burnin=100000, n.sample=300000, thin=100
)

```

`hex.m` is a 914 $\times$ 914 neighbourhood matrix with its elements assuming the value 1 if the two corresponding hexagon cells share a boundary and 0 otherwise. The formula is run to generate three independent Markov chains, which are used for the inferential analysis. Each times it runs the formula will generate 300,000 MCMC samples (`n.sample`), 100,000 of which are then discarded (the so-called `burnin` period), while the remaining 200,000 are `thin`ned by 100 to reduce their temporal autocorrelation^[Lee, Duncan. 2013. "Carbayes: An R Package for Bayesian Spatial Modeling with Conditional Autoregressive Priors." Journal of Statistical Software 55 (13): 1–24. https://doi.org/10.18637/jss.v055.i13.].  

## Variables

```{r}

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-prepare",
                    object = 'vars_1.df',
                    path = knitr_cache_path),
  file='NUL')

vars_1.df %>%
  dplyr::mutate(`Meetup RSVPs Period 1 2007-12 (%)` = rsvp_sum_2007_2012_perc*100,
                `Meetup RSVPs Period 1 2007-12 (sqrt)` = sqrt(rsvp_sum_2007_2012_perc),
                `Meetup RSVPs Period 2 2012-13 (%)` = rsvp_sum_2012_2013_perc*100,
                `Meetup RSVPs Period 2 2012-13 (sqrt)` = sqrt(rsvp_sum_2012_2013_perc),
                `Meetup RSVPs Period 2 (sqrt) - Period 1 (sqrt)` = rsvp_sum_diff*100,
                `hh with Internet 2007-12` = dd_2007_2012_perc*100,
                `hh with Internet 2013` = dd_2013_perc*100,
                `M5S vote 2013 (%)` = M5S_2013_perc*100,
                `turnout 2006-08 (mean %)` = turnout_2006_08*100,
                `staff not-for-profit ass. (%)` = NP_949_perc*100,
                `staff not-for-profit ass. (sqrt)` = sqrt(NP_949_perc),
                `pop. with degree (%)` = P47_perc*100,
                `pop. unemployed (%)` = P62_perc*100,
                `pop. homemaker (%)` = P130_perc*100,
                `pop. over 65 (%)` = over65_perc*100) %>%
  dplyr::select(`Meetup RSVPs Period 1 2007-12 (%)`:`pop. over 65 (%)`) %>%
  tidyr::pivot_longer(cols = `Meetup RSVPs Period 1 2007-12 (%)`:`pop. over 65 (%)`) %>%
  dplyr::group_by(name) %>% 
  summarize_all(funs(min, max, mean, median, sd), na.rm = T) %>%
  kbl(booktabs = T, caption = "Spatial analysis regression models: Continuous variables", digits=4) %>% 
  kable_styling(latex_options = c("striped", "HOLD_position"), font_size = 7) %>%
  column_spec(1, bold = T, width = "25em")

```

```{r fig.cap = "Distribution of continuous for used in the spatial analysis", fig.width = 20, fig.height = 20}

tmp.df <- 
  vars_1.df %>%
  dplyr::mutate(`Meetup RSVPs Period 1 2007-12 (%)` = rsvp_sum_2007_2012_perc*100,
                `Meetup RSVPs Period 1 2007-12 (sqrt)` = sqrt(rsvp_sum_2007_2012_perc),
                `Meetup RSVPs Period 2 2012-13 (%)` = rsvp_sum_2012_2013_perc*100,
                `Meetup RSVPs Period 2 2012-13 (sqrt)` = sqrt(rsvp_sum_2012_2013_perc),
                `Meetup RSVPs Period 2 (sqrt) - Period 1 (sqrt)` = rsvp_sum_diff*100,
                `hh with Internet 2007-12` = dd_2007_2012_perc*100,
                `hh with Internet 2013` = dd_2013_perc*100,
                `M5S vote 2013 (%)` = M5S_2013_perc*100,
                `turnout 2006-08 (mean %)` = turnout_2006_08*100,
                `staff not-for-profit ass. (%)` = NP_949_perc*100,
                `staff not-for-profit ass. (sqrt)` = sqrt(NP_949_perc),
                `pop. with degree (%)` = P47_perc*100,
                `pop. unemployed (%)` = P62_perc*100,
                `pop. homemaker (%)` = P130_perc*100,
                `pop. over 65 (%)` = over65_perc*100) %>%
  dplyr::select(hex, `Meetup RSVPs Period 1 2007-12 (%)`:`pop. over 65 (%)`)

tmp.gg <- list()

for (i in 2:ncol(tmp.df)) {
  
  tmp.sf <- 
    italy_hex_cropped.sf
  
  tmp.sf$this_var <- 
    tmp.df[[i]][match(tmp.sf$index, tmp.df$hex)]
  
  tmp.gg[[i]] <- 
    cowplot::plot_grid(
      ggplot(tmp.sf, aes(fill = this_var)) +
        geom_sf(colour = NA) + guides(fill = 'none') +
        scale_fill_distiller(palette = "YlGnBu") +
        ggplot2::theme_bw() + 
        ggplot2::theme(axis.text = ggplot2::element_blank(), 
                       axis.ticks = ggplot2::element_blank(), 
                       axis.title = ggplot2::element_blank(), 
                       legend.key = ggplot2::element_blank(), 
                       panel.background = ggplot2::element_rect(fill = "white", 
                                                                colour = NA), 
                       panel.border =  ggplot2::element_blank(), 
                       panel.grid = ggplot2::element_blank()),
      ggplot(tmp.sf, aes(this_var)) +
        geom_density() +
        labs(x = names(tmp.df)[i]),
      ncol = 1, 
      rel_heights = c(1, .3)
    )
  
}

cowplot::plot_grid(plotlist = tmp.gg[c(2:5, 7:8, 9:14, 6)],
                   ncol = 5)

```

```{r}

vars_1.df %>% 
  dplyr::mutate(`pop. density (ordinal)` = P1_cut_6) %>%
  dplyr::select(`pop. density (ordinal)`) %>%
  summary() %>%
  t() %>%
  kbl(booktabs = T, caption = "Categorical variable", digits=4) %>% 
  kable_styling(latex_options = c("striped", "HOLD_position"), font_size = 7) %>%
  column_spec(1, bold = T, width = "15em")

``` 

## Model summaries

```{r}

spatial_models_variable_map <- 
  c(rsvp_sum_2007_2012_perc = "Meetup RSVPs Period 1 2007-12",
    rsvp_sum_2012_2013_perc = "Meetup RSVPs Period 2 2012-13",
    rsvp_sum_diff = "Meetup RSVPs Period 2 - Period 1",
    dd_2007_2012_perc = "hh with Internet 2007-12",
    dd_2013_perc = "hh with Internet 2013",
    M5S_2013_perc = "M5S vote 2013",
    turnout_2006_08 = "turnout 2006-08",
    NP_949_perc = "staff not-for-profit ass.",
    P47_perc = "pop. with degree",
    P62_perc = "pop. unemployed",
    P130_perc = "pop. homemaker",
    over65_perc = "pop. over 65",
    P1_cut_6 = "pop. density")

```

### Predicting meetup RSVPs in 2007-2012: Model 1

```{r}

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-1",
                    object = 'var_1_MCMC.chain_1',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-1",
                    object = 'var_1_MCMC.chain_2',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-1",
                    object = 'var_1_MCMC.chain_3',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-1",
                    object = 'beta_var_1.samples_2007_2012',
                    path = knitr_cache_path),
  file='NUL')

```

```{r}

var_1_MCMC.chain_1$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                           "turnout 2006-08",
                           "staff not-for-profit ass. (sqrt)",
                           "hh with Internet",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_1_MCMC.chain_1$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting meetup RSVPs in 2007-2012: Model 1, Chain 1") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

```{r}

var_1_MCMC.chain_2$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                           "turnout 2006-08",
                           "staff not-for-profit ass. (sqrt)",
                           "hh with Internet",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_1_MCMC.chain_2$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting meetup RSVPs in 2007-2012: Model 1, Chain 2") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

```{r}

var_1_MCMC.chain_3$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                           "turnout 2006-08",
                           "staff not-for-profit ass. (sqrt)",
                           "hh with Internet",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_1_MCMC.chain_3$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting meetup RSVPs in 2007-2012: Model 1, Chain 3") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

```{r fig.pos = "H", fig.cap = "Predicting meetup RSVPs in 2007-2012: Model 1, Convergence of the three Markov chains", fig.align = 'left', fig.width = 10, fig.height = 6}

plot(beta_var_1.samples_2007_2012[ ,2:4])

```

```{r}

coda::gelman.diag(beta_var_1.samples_2007_2012)$psrf %>%
  kbl(booktabs = T, 
      caption = 
        sprintf("Predicting meetup RSVPs in 2007-2012: Model 1, Potential scale reduction factor (multivariate PSRF: %s)",
                round(coda::gelman.diag(beta_var_1.samples_2007_2012)$mpsrf, 6))) %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

### Predicting meetup RSVPs in 2007-2012: Model 2

```{r}

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-1",
                    object = 'var_2_MCMC.chain_1',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-1",
                    object = 'var_2_MCMC.chain_2',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-1",
                    object = 'var_2_MCMC.chain_3',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-1",
                    object = 'beta_var_2.samples_2007_2012',
                    path = knitr_cache_path),
  file='NUL')

```


```{r}

var_2_MCMC.chain_1$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                           "hh with Internet",
                           "pop. with degree",
                           "pop. unemployed",
                           "pop. homemaker",
                           "pop. over 65",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_2_MCMC.chain_1$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting meetup RSVPs in 2007-2012: Model 2, Chain 1") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

```{r}

var_2_MCMC.chain_2$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                          "hh with Internet",
                           "pop. with degree",
                           "pop. unemployed",
                           "pop. homemaker",
                           "pop. over 65",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_2_MCMC.chain_2$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting meetup RSVPs in 2007-2012: Model 2, Chain 2") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

```{r}

var_2_MCMC.chain_3$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                          "hh with Internet",
                           "pop. with degree",
                           "pop. unemployed",
                           "pop. homemaker",
                           "pop. over 65",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_2_MCMC.chain_3$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting meetup RSVPs in 2007-2012: Model 2, Chain 3") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

```{r fig.pos = "H", fig.cap = "Predicting meetup RSVPs in 2007-2012: Model 2, Convergence of the three Markov chains", fig.align = 'left', fig.width = 10, fig.height = 6}

plot(beta_var_2.samples_2007_2012[ ,2:4])

```

```{r}

coda::gelman.diag(beta_var_2.samples_2007_2012)$psrf %>%
  kbl(booktabs = T, 
      caption = 
        sprintf("Predicting meetup RSVPs in 2007-2012: Model 2, Potential scale reduction factor (multivariate PSRF: %s)",
                round(coda::gelman.diag(beta_var_2.samples_2007_2012)$mpsrf, 6))) %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

### Predicting meetup RSVPs in 2012-2013: Model 1

```{r}

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-2",
                    object = 'var_1_MCMC.chain_1',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-2",
                    object = 'var_1_MCMC.chain_2',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-2",
                    object = 'var_1_MCMC.chain_3',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-2",
                    object = 'beta_var_1.samples_2012_2013',
                    path = knitr_cache_path),
  file='NUL')

```

```{r}

var_1_MCMC.chain_1$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                           "turnout 2006-08",
                           "staff not-for-profit ass. (sqrt)",
                           "hh with Internet",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_1_MCMC.chain_1$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting meetup RSVPs in 2012-2013: Model 1, Chain 1") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

```{r}

var_1_MCMC.chain_2$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                           "turnout 2006-08",
                           "staff not-for-profit ass. (sqrt)",
                           "hh with Internet",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_1_MCMC.chain_2$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting meetup RSVPs in 2012-2013: Model 1, Chain 2") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

```{r}

var_1_MCMC.chain_3$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                           "turnout 2006-08",
                           "staff not-for-profit ass. (sqrt)",
                           "hh with Internet",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_1_MCMC.chain_3$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting meetup RSVPs in 2012-2013: Model 1, Chain 3") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

```{r fig.pos = "H", fig.cap = "Predicting meetup RSVPs in 2012-2013: Model 1, Convergence of the three Markov chains", fig.align = 'left', fig.width = 10, fig.height = 6}

plot(beta_var_1.samples_2012_2013[ ,2:4])

```

```{r}

coda::gelman.diag(beta_var_1.samples_2012_2013)$psrf %>%
  kbl(booktabs = T, 
      caption = 
        sprintf("Predicting meetup RSVPs in 2012-2013: Model 1, Potential scale reduction factor (multivariate PSRF: %s)",
                round(coda::gelman.diag(beta_var_1.samples_2012_2013)$mpsrf, 6))) %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

### Predicting meetup RSVPs in 2012-2013: Model 2

```{r}

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-2",
                    object = 'var_2_MCMC.chain_1',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-2",
                    object = 'var_2_MCMC.chain_2',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-2",
                    object = 'var_2_MCMC.chain_3',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-2",
                    object = 'beta_var_2.samples_2012_2013',
                    path = knitr_cache_path),
  file='NUL')

```

```{r}

var_2_MCMC.chain_1$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                           "hh with Internet",
                           "pop. with degree",
                           "pop. unemployed",
                           "pop. homemaker",
                           "pop. over 65",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_2_MCMC.chain_1$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting meetup RSVPs in 2012-2013: Model 2, Chain 1") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

```{r}

var_2_MCMC.chain_2$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                           "hh with Internet",
                           "pop. with degree",
                           "pop. unemployed",
                           "pop. homemaker",
                           "pop. over 65",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_2_MCMC.chain_2$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting meetup RSVPs in 2012-2013: Model 2, Chain 2") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

```{r}

var_2_MCMC.chain_3$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                           "hh with Internet",
                           "pop. with degree",
                           "pop. unemployed",
                           "pop. homemaker",
                           "pop. over 65",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_2_MCMC.chain_3$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting meetup RSVPs in 2012-2013: Model 2, Chain 3") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

```{r fig.pos = "H", fig.cap = "Predicting meetup RSVPs in 2012-2013: Model 2, Convergence of the three Markov chains", fig.align = 'left', fig.width = 10, fig.height = 6}

plot(beta_var_2.samples_2012_2013[ ,2:4])

```

```{r}

coda::gelman.diag(beta_var_2.samples_2012_2013)$psrf %>%
  kbl(booktabs = T, 
      caption = 
        sprintf("Predicting meetup RSVPs in 2012-2013: Model 2, Potential scale reduction factor (multivariate PSRF: %s)",
                round(coda::gelman.diag(beta_var_2.samples_2012_2013)$mpsrf, 6))) %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

### Predicting Meetup RSVPs Period 2 - Period 1: Model 1

```{r}

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-5",
                    object = 'var_6_MCMC.chain_1',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-5",
                    object = 'var_6_MCMC.chain_2',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-5",
                    object = 'var_6_MCMC.chain_3',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-5",
                    object = 'beta_var_6.samples_difff',
                    path = knitr_cache_path),
  file='NUL')

```

```{r}

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-2",
                    object = 'var_6_MCMC.chain_1',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-2",
                    object = 'var_6_MCMC.chain_2',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-2",
                    object = 'var_6_MCMC.chain_3',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-2",
                    object = 'beta_var_6.samples_diff',
                    path = knitr_cache_path),
  file='NUL')

```

```{r}

var_6_MCMC.chain_1$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                           "turnout 2006-08",
                           "staff not-for-profit ass. (sqrt)",
                           "hh with Internet",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_6_MCMC.chain_1$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting Meetup RSVPs Period 2 - Period 1 : Model 1, Chain 1") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

```{r}

var_6_MCMC.chain_2$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                           "turnout 2006-08",
                           "staff not-for-profit ass. (sqrt)",
                           "hh with Internet",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_6_MCMC.chain_2$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting Meetup RSVPs Period 2 - Period 1: Model 1, Chain 2") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

```{r}

var_6_MCMC.chain_3$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                           "turnout 2006-08",
                           "staff not-for-profit ass. (sqrt)",
                           "hh with Internet",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_6_MCMC.chain_3$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting Meetup RSVPs Period 2 - Period 1: Model 1, Chain 3") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

```{r fig.pos = "H", fig.cap = "Predicting meetup RSVPs in 2012-2013: Model 1, Convergence of the three Markov chains", fig.align = 'left', fig.width = 10, fig.height = 6}

plot(beta_var_6.samples_diff[ ,2:4])

```

```{r}

coda::gelman.diag(beta_var_6.samples_diff)$psrf %>%
  kbl(booktabs = T, 
      caption = 
        sprintf("Predicting Meetup RSVPs Period 2 - Period 1: Model 1, Potential scale reduction factor (multivariate PSRF: %s)",
                round(coda::gelman.diag(beta_var_6.samples_diff)$mpsrf, 6))) %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

### Predicting Meetup RSVPs Period 2 - Period 1: Model 2

```{r}

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-5",
                    object = 'var_7_MCMC.chain_1',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-5",
                    object = 'var_7_MCMC.chain_2',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-5",
                    object = 'var_7_MCMC.chain_3',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-5",
                    object = 'beta_var_7.samples_difff',
                    path = knitr_cache_path),
  file='NUL')

```

```{r}

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-2",
                    object = 'var_7_MCMC.chain_1',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-2",
                    object = 'var_7_MCMC.chain_2',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-2",
                    object = 'var_7_MCMC.chain_3',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-2",
                    object = 'beta_var_7.samples_diff',
                    path = knitr_cache_path),
  file='NUL')

```

```{r}

var_7_MCMC.chain_1$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                           "hh with Internet",
                           "pop. with degree",
                           "pop. unemployed",
                           "pop. homemaker",
                           "pop. over 65",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_7_MCMC.chain_1$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting Meetup RSVPs Period 2 - Period 1: Model 2, Chain 1") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

```{r}

var_7_MCMC.chain_2$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                           "hh with Internet",
                           "pop. with degree",
                           "pop. unemployed",
                           "pop. homemaker",
                           "pop. over 65",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_7_MCMC.chain_2$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting Meetup RSVPs Period 2 - Period 1: Model 2, Chain 2") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

```{r}

var_7_MCMC.chain_3$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                           "hh with Internet",
                           "pop. with degree",
                           "pop. unemployed",
                           "pop. homemaker",
                           "pop. over 65",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_7_MCMC.chain_3$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting Meetup RSVPs Period 2 - Period 1: Model 2, Chain 3") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

```{r fig.pos = "H", fig.cap = "Predicting meetup RSVPs in 2012-2013: Model 2, Convergence of the three Markov chains", fig.align = 'left', fig.width = 10, fig.height = 6}

plot(beta_var_7.samples_diff[ ,2:4])

```

```{r}

coda::gelman.diag(beta_var_7.samples_diff)$psrf %>%
  kbl(booktabs = T, 
      caption = 
        sprintf("Predicting Meetup RSVPs Period 2 - Period 1: Model 2, Potential scale reduction factor (multivariate PSRF: %s)",
                round(coda::gelman.diag(beta_var_7.samples_diff)$mpsrf, 6))) %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

### Predicting meetup RSVPs with 7-day rolling windows

```{r}

load("data/out/meetup_NP_949_rolling_posts.df.RData")

tmp.df <- 
  rolling_posts.df %>%
  tibble::rownames_to_column("var") 

tmp.df$var <-
  gsub("turnout_2006_08", "turnout 2006-08", tmp.df$var)

tmp.df$var <-
  gsub("sqrt\\(NP_949_perc\\)", "staff not-for-profit ass. (sqrt)", tmp.df$var)

tmp.df$var <-
  gsub("dd_perc", "hh with Internet", tmp.df$var)

tmp.df$var <-
  gsub("P1_cut_6\\.L", "pop. density (ordinal, linear)", tmp.df$var)

tmp.df$var <-
  gsub("P1_cut_6\\.Q", "pop. density (ordinal, quadratic)", tmp.df$var)

tmp.df$var <-
  gsub("P1_cut_6\\.C", "pop. density (ordinal, cubic)", tmp.df$var)

tmp.df$var <-
  gsub("P1_cut_6\\^4", "pop. density (ordinal, quartic)", tmp.df$var)

tmp.df %>%  
  dplyr::mutate(var = gsub("\\.\\.\\.(.*)", "", var),
                week = as.Date(week)) %>%
  dplyr::group_by(var) %>%
  dplyr::summarise(From = min(week),
                   To = max(week),
                   `Mean coefficient` = mean(V1),
                   `Mean 2.5%` = mean(`2.5%`),
                   `Mean 97.5%` = mean(`97.5%`),
                   Num.Obs. = n(),
                   `Num.Sign. >0` = sum(`2.5%` > 0),
                   `Perc.Sign. >0` = sum(`2.5%` > 0) / n() * 100,
                   `Num.Sign. <0` = sum(`97.5%` < 0),
                   `Perc.Sign. <0` = sum(`97.5%` < 0) / n() * 100) %>%
  
  dplyr::bind_rows(
    tmp.df %>%  
      dplyr::filter(week >= break_dates[1] & week < break_dates[2]) %>%
      dplyr::mutate(var = gsub("\\.\\.\\.(.*)", "", var),
                    week = as.Date(week)) %>%
      dplyr::group_by(var) %>%
      dplyr::summarise(From = min(week),
                       To = max(week),
                       `Mean coefficient` = mean(V1),
                       `Mean 2.5%` = mean(`2.5%`),
                       `Mean 97.5%` = mean(`97.5%`),
                       Num.Obs. = n(),
                       `Num.Sign. >0` = sum(`2.5%` > 0),
                       `Perc.Sign. >0` = sum(`2.5%` > 0) / n() * 100,
                       `Num.Sign. <0` = sum(`97.5%` < 0),
                       `Perc.Sign. <0` = sum(`97.5%` < 0) / n() * 100)
  ) %>%
  
  dplyr::bind_rows(
    tmp.df %>%  
      dplyr::filter(week >= break_dates[2] & week < break_dates[3]) %>%
      dplyr::mutate(var = gsub("\\.\\.\\.(.*)", "", var),
                    week = as.Date(week)) %>%
      dplyr::group_by(var) %>%
      dplyr::summarise(From = min(week),
                       To = max(week),
                       `Mean coefficient` = mean(V1),
                       `Mean 2.5%` = mean(`2.5%`),
                       `Mean 97.5%` = mean(`97.5%`),
                       Num.Obs. = n(),
                       `Num.Sign. >0` = sum(`2.5%` > 0),
                       `Perc.Sign. >0` = sum(`2.5%` > 0) / n() * 100,
                       `Num.Sign. <0` = sum(`97.5%` < 0),
                       `Perc.Sign. <0` = sum(`97.5%` < 0) / n() * 100)
  ) %>%
  
  kbl(booktabs = T, 
      caption = "Means of coefficients and confidence interval boundaries predicting Meetup RVSPs (weekly count)", digits = 2) %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7) %>%
  kableExtra::column_spec(1, width = "2.5cm") %>%
  kableExtra::column_spec(4:11, width = ".8cm") %>%
  kableExtra::row_spec(8, hline_after = T) %>%
  kableExtra::row_spec(16, hline_after = T)

```


### Predicting M5S vote in 2013 with Meetup RSVPs Period 1 2007-12

```{r}

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-3",
                    object = 'var_3_MCMC.chain_1',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-3",
                    object = 'var_3_MCMC.chain_2',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-3",
                    object = 'var_3_MCMC.chain_3',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-3",
                    object = 'beta_var_3.samples',
                    path = knitr_cache_path),
  file='NUL')

```

```{r}

var_3_MCMC.chain_1$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                           "Meetup RSVPs Period 1 2007-12 (sqrt)",
                           "hh with Internet",
                           "pop. with degree",
                           "pop. unemployed",
                           "pop. homemaker",
                           "pop. over 65",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_3_MCMC.chain_1$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting M5S vote in 2013 with Meetup RSVPs Period 1 2007-12, Chain 1") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

```{r}

var_3_MCMC.chain_2$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                           "Meetup RSVPs Period 1 2007-12 (sqrt)",
                           "hh with Internet",
                           "pop. with degree",
                           "pop. unemployed",
                           "pop. homemaker",
                           "pop. over 65",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_3_MCMC.chain_2$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting M5S vote in 2013 with Meetup RSVPs Period 1 2007-12, Chain 2") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

```{r}

var_3_MCMC.chain_3$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                           "Meetup RSVPs Period 1 2007-12 (sqrt)",
                           "hh with Internet",
                           "pop. with degree",
                           "pop. unemployed",
                           "pop. homemaker",
                           "pop. over 65",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_3_MCMC.chain_3$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting M5S vote in 2013 with Meetup RSVPs Period 1 2007-12, Chain 3") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

```{r fig.pos = "H", fig.cap = "Predicting M5S vote in 2013 with Meetup RSVPs Period 1 2007-12, Convergence of the three Markov chains", fig.align = 'left', fig.width = 10, fig.height = 6}

plot(beta_var_3.samples[ ,2:4])

```

```{r}

gelman.diag(beta_var_3.samples)$psrf %>%
  kbl(booktabs = T, 
      caption = 
        sprintf("Predicting M5S vote in 2013 with Meetup RSVPs Period 1 2007-12, Potential scale reduction factor (multivariate PSRF: %s)",
                round(coda::gelman.diag(beta_var_3.samples)$mpsrf, 6))) %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

### Predicting M5S vote in 2013 with Meetup RSVPs Period 2 2012-13

```{r}

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-3",
                    object = 'var_4_MCMC.chain_1',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-3",
                    object = 'var_4_MCMC.chain_2',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-3",
                    object = 'var_4_MCMC.chain_3',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-3",
                    object = 'beta_var_4.samples',
                    path = knitr_cache_path),
  file='NUL')

```

```{r}

var_4_MCMC.chain_1$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                           "Meetup RSVPs Period 2 2012-13 (sqrt)",
                           "hh with Internet",
                           "pop. with degree",
                           "pop. unemployed",
                           "pop. homemaker",
                           "pop. over 65",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_4_MCMC.chain_1$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting M5S vote in 2013 with Meetup RSVPs Period 2 2012-13, Chain 1") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

```{r}

var_4_MCMC.chain_2$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                           "Meetup RSVPs Period 2 2012-13 (sqrt)",
                           "hh with Internet",
                           "pop. with degree",
                           "pop. unemployed",
                           "pop. homemaker",
                           "pop. over 65",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_4_MCMC.chain_2$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting M5S vote in 2013 with Meetup RSVPs Period 2 2012-13, Chain 2") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```


```{r}

var_4_MCMC.chain_3$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                           "Meetup RSVPs Period 2 2012-13 (sqrt)",
                           "hh with Internet",
                           "pop. with degree",
                           "pop. unemployed",
                           "pop. homemaker",
                           "pop. over 65",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_4_MCMC.chain_3$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting M5S vote in 2013 with Meetup RSVPs Period 2 2012-13, Chain 3") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

```{r fig.pos = "H", fig.cap = "Predicting M5S vote in 2013 with Meetup RSVPs Period 2 2012-13, Convergence of the three Markov chains", fig.align = 'left', fig.width = 10, fig.height = 6}

plot(beta_var_4.samples[ ,2:4])

```

```{r}

gelman.diag(beta_var_4.samples)$psrf %>%
  kbl(booktabs = T, 
      caption = 
        sprintf("Predicting M5S vote in 2013 with Meetup RSVPs Period 2 2012-13, Potential scale reduction factor (multivariate PSRF: %s)",
                round(coda::gelman.diag(beta_var_4.samples)$mpsrf, 6))) %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```

### Predicting M5S vote in 2013 with Meetup RSVPs Period 2 - Period 1 

```{r}

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-3",
                    object = 'var_5_MCMC.chain_1',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-3",
                    object = 'var_5_MCMC.chain_2',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-3",
                    object = 'var_5_MCMC.chain_3',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "meetup-spatial-analysis-model-3",
                    object = 'beta_var_5.samples',
                    path = knitr_cache_path),
  file='NUL')

```

```{r}

var_5_MCMC.chain_1$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                           "Meetup RSVPs Period 2 (sqrt) - Period 1 (sqrt)",
                           "hh with Internet",
                           "pop. with degree",
                           "pop. unemployed",
                           "pop. homemaker",
                           "pop. over 65",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_5_MCMC.chain_1$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting M5S vote in 2013 with Meetup RSVPs Period 2 - Period 1, Chain 1") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7) %>%
  column_spec(1, width = "4 cm")

```

```{r}

var_5_MCMC.chain_2$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                           "Meetup RSVPs Period 2 (sqrt) - Period 1 (sqrt)",
                           "hh with Internet",
                           "pop. with degree",
                           "pop. unemployed",
                           "pop. homemaker",
                           "pop. over 65",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_5_MCMC.chain_2$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting M5S vote in 2013 with Meetup RSVPs Period 2 - Period 1, Chain 2") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7) %>%
  column_spec(1, width = "4 cm")

```


```{r}

var_5_MCMC.chain_3$summary.results %>%
  magrittr::set_rownames(c("(Intercept)",
                           "Meetup RSVPs Period 2 (sqrt) - Period 1 (sqrt)",
                           "hh with Internet",
                           "pop. with degree",
                           "pop. unemployed",
                           "pop. homemaker",
                           "pop. over 65",
                           "pop. density (ordinal, linear)",
                           "pop. density (ordinal, quadratic)",
                           "pop. density (ordinal, cubic)",
                           "pop. density (ordinal, quartic)",
                           "nu2",
                           "tau2",
                           "rho" )) %>%
  as.data.frame() %>%
  dplyr::mutate(`Num.Obs.` = length(var_5_MCMC.chain_3$fitted.values)) %>%
  kbl(booktabs = T, 
      caption = "Predicting M5S vote in 2013 with Meetup RSVPs Period 2 - Period 1, Chain 3") %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7) %>%
  column_spec(1, width = "4 cm")

```

```{r fig.pos = "H", fig.cap = "Predicting M5S vote in 2013 with Meetup RSVPs Period 2 - Period 1, Convergence of the three Markov chains", fig.align = 'left', fig.width = 10, fig.height = 6}

plot(beta_var_5.samples[ ,2:4])

```

```{r}

gelman.diag(beta_var_5.samples)$psrf %>%
  kbl(booktabs = T, 
      caption = 
        sprintf("Predicting M5S vote in 2013 with Meetup RSVPs Period 2 - Period 1, Potential scale reduction factor (multivariate PSRF: %s)",
                round(coda::gelman.diag(beta_var_5.samples)$mpsrf, 6))) %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 7)

```


# Survey responses: confidence intervals and regression models

I specified three models to conduct the analysis at the individual level based on survey responses collected during the 2013 election study:^[Associazione Itanes. 2013. “Italian National Election Studies.” 2013. http://www.itanes.org/dati/]

1. A model to predict political discontent with external political efficacy among the independent variables.

2. A model to predict internal political efficacy with the difference in meetup activity between the two periods among the independent variables. 

3. A model to predict the M5S vote with an interaction term including political discontent and internal political efficacy. 

All independent continuous variables have been scaled to facilitate the visual interpretation of the coefficients. 

## Variables

```{r}

individual_variable_mapping_cont <- 
  c("rsvp_sum_2007_2012_perc" = "Meetup RSVPs Period 1 (per 10,000 ppl)",
    "rsvp_sum_2012_2013_perc" = "Meetup RSVPs Period 2 (per 10,000 ppl)",
    "rsvp_sum_diff" = "Meetup RSVPs Period 2 (sqrt) - Period 1 (sqrt)",
    "age" = "age",
    "education" = "education",
    "municipality_size" = "municipality size",
    "pol_efficacy_int" = "freq. political talks",
    "pol_efficacy_ext" = "ext. political efficacy",
    "NOSAY" = "NOSAY",
    "NOCARE" = "NOCARE",
    "DISTANCE" = "DISTANCE",
    "COMPLEXITY" = "COMPLEXITY",
    "pol_interest" = "political interest",
    "pol_part" = "political participation",
    "discontent" = "discontent")

individual_variable_mapping_cat <- 
  c("voted_M5S_2013" = "voted M5S in 2013",
    "sex" = "sex",
    "social_capital_bin" = "member of association")

```

```{r}

individual_variable_questions <- 
  c("age" = 
      "D1. Could you please tell me your age? [numerical]",
    "education" = 
      "STU. What’s your level of education? [numerical]",
    "municipality_size" = 
      "[inferred] size of the municipality of residence",
    "pol_efficacy_int" = 
      "D28. How often have you discussed political matters during the last two months? [numerical]",
    "NOSAY" = 
      "D38_1. People like me don’t have any say in what the government does (what degree they are true, in your opinion?) [numerical]",
    "NOCARE" = 
      "D38_4. Parties are only interested in people's votes, but not in their opinions  (what degree they are true, in your opinion?) [numerical]",
    "DISTANCE" = 
      "D38_3. Usually, people we elect to the Parliament quickly lose touch with the people  (what degree they are true, in your opinion?) [numerical]",
    "COMPLEXITY" = 
      "D38_2. Sometimes politics and government seem so complicated that a person like me can’t really understand what’s going on (what degree they are true, in your opinion?) [numerical]",
    "pol_efficacy_ext" = 
      "NOSAY + NOCARE + DISTANCE + COMPLEXITY [numerical]",
    "pol_interest" = 
      "D27: In general, are you interested in politics? [numerical]",
    "pol_part" = 
      "D78. I am now going to read a list of different associations. For each of them, please tell me whether you are a member and/or participate in their activities. 1. Voluntary associations + 2. Religious associations; 3. Sports associations + 4. Cultural associations + 5. Environmental associations + 6. Third-world and human rights associations + 7. Worker unions/trade associations + 8. Movements or political parties [numerical]",
    "discontent" = 
      "D15. I am going to read you a list of institutions. Can you tell me how much confidence you have in them: none at all, not very much, a fair amount or a great deal? D15_1. The Parliament + D15_2. The Political Parties [numerical]",
    "voted_M5S_2013" = 
      "D90. Can you tell me which party you voted for on the Chamber of Deputies ballot? [binary]",
    "sex" = 
      "[not asked, inferred, binary]",
    "social_capital_bin" = 
      "D78. I am now going to read a list of different associations. For each of them, please tell me whether you are a member and/or participate in their activities. 1. Voluntary associations OR 5. Environmental associations OR 6. Third-world and human rights associationss OR 8. Movements or political parties [binary]")

data.frame(label = c(names(individual_variable_mapping_cont),
                     names(individual_variable_mapping_cat)),
           variable = c(individual_variable_mapping_cont,
                        individual_variable_mapping_cat)) %>%
  dplyr::inner_join(data.frame(label = names(individual_variable_questions),
                               question = individual_variable_questions),
                    by = 'label') %>%
  dplyr::select(variable, question) %>%
  kbl(booktabs = T, caption = "Survey questions and coding of variables") %>% 
  kable_styling(latex_options = c("striped", "HOLD_position"), font_size = 7) %>%
  column_spec(1, bold = T, width = "3 cm") %>%
  column_spec(2, bold = F, width = "15 cm")

```

```{r}
capture.output(
  knitr::load_cache(label = "itanes-individual-level",
                    object = 'discontent.lm_1',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "itanes-individual-level",
                    object = 'discontent.lm_2',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "itanes-individual-level",
                    object = 'pol_efficacy_int_1.lmer',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "itanes-individual-level",
                    object = 'pol_efficacy_int_1a.lmer',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "itanes-individual-level",
                    object = 'pol_efficacy_int_1b.lmer',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "itanes-individual-level",
                    object = 'pol_efficacy_int_2.lmer',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "itanes-individual-level",
                    object = 'pol_efficacy_int_2a.lmer',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "itanes-individual-level",
                    object = 'pol_efficacy_int_2b.lmer',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "itanes-individual-level",
                    object = 'pol_efficacy_int_3.lm',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "itanes-individual-level",
                    object = 'pol_efficacy_int_4.lm',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "itanes-individual-level",
                    object = 'voted_m5s_2013.glm_1',
                    path = knitr_cache_path),
  file='NUL')

capture.output(
  knitr::load_cache(label = "itanes-individual-level",
                    object = 'voted_m5s_2013.glm_2',
                    path = knitr_cache_path),
  file='NUL')

```


## Model summaries

```{r}

individual_variable_mapping_regressions <- 
  c("rsvp_sum_2007_2012_perc" = "Meetup RSVPs Period 1",
    "rsvp_sum_2012_2013_perc" = "Meetup RSVPs Period 2",
    "scale(sqrt(rsvp_sum_2007_2012_perc))" = "Meetup RSVPs Period 1",
    "scale(sqrt(rsvp_sum_2012_2013_perc))" = "Meetup RSVPs Period 2",
    "rsvp_sum_diff" = "Meetup RSVPs Period 2 (sqrt) - Period 1 (sqrt)",
    "age" = "age",
    "education" = "education",
    "municipality_size" = "municipality size",
    "pol_efficacy_int" = "freq. political talks",
    "pol_efficacy_ext" = "ext. political efficacy",
    "NOSAY" = "NOSAY",
    "NOCARE" = "NOCARE",
    "DISTANCE" = "DISTANCE",
    "COMPLEXITY" = "COMPLEXITY",
    "pol_interest" = "political interest",
    "pol_part" = "political participation",
    "discontent" = "discontent",
    "voted_M5S_2013TRUE" = "voted M5S in 2013",
    "sexFemale" = "sex: Female",
    "social_capital_binTRUE" = "member of association",
    "pol_efficacy_int:discontent" = "freq. political talks  X discontent")

```

Notes to Table \@ref(tab:itanes-tab-lm): 

* LM, mod. 1 and mod. 2 are not described in the paper. 

* In LM, mod. 2 the external political efficacy index is replaced by its four components.

```{r itanes-tab-lm}

modelsummary(list("mod. 1" = discontent.lm_1, 
                  "mod. 2" = discontent.lm_2), 
             output = 'kableExtra',
             coef_map = individual_variable_mapping_regressions, stars = TRUE, 
             gof_omit = 'IC|Log|Adj',
             caption = "Linear regression models predicting political discontent using survey responses")  %>%
  kable_styling(font_size = 7) %>%
  add_header_above(c(" " = 1, "LM" = 2))

```

Notes to Table \@ref(tab:itanes-tab-lme): 

* LME, Mod. 1 (\*) is described in the paper.

* LME, mod. 1-6 are linear mixed-effects models with fixed mean random intercepts and municipality as the grouping factor `g` (with formula `(1 | g)`).

* LM, mod. 1-2 are simple linear models (with no hierarchical structure).

* LME, mod. 1-3 and LM mod. 1 exclude from the analysis respondents from the metropolitan municipalities of Rome, Milan and Naples.

```{r itanes-tab-lme}

modelsummary(list("mod. 1 (*)" = pol_efficacy_int_1.lmer,
                  "mod. 2" = pol_efficacy_int_1a.lmer,
                  "mod. 3" = pol_efficacy_int_1b.lmer,
                  "mod. 4" = pol_efficacy_int_2.lmer,
                  "mod. 5" = pol_efficacy_int_2a.lmer,
                  "mod. 6" = pol_efficacy_int_2b.lmer,
                  "mod. 1" = pol_efficacy_int_3.lm,
                  "mod. 2" = pol_efficacy_int_4.lm), 
             output = 'kableExtra',
             coef_map = individual_variable_mapping_regressions, stars = TRUE, 
             gof_omit = 'IC|Log|Adj',
             caption = "Linear regression models predicting freq. of political talks using survey responses") %>%
  add_header_above(c(" " = 1, "LME" = 6, "LM" = 2)) %>%
  kable_styling(font_size = 7) %>%
  column_spec(1, bold = F, width = "3.5cm") 

```

Notes to Table \@ref(tab:itanes-tab-bin): 

* BIN mod. 2 (\*) is described in the paper.

```{r itanes-tab-bin}

modelsummary(list("mod. 1" = voted_m5s_2013.glm_1,
                  "mod. 2 (*)" = voted_m5s_2013.glm_2), 
             output = 'kableExtra',
             coef_map = individual_variable_mapping_regressions, stars = TRUE, 
             gof_omit = 'IC|Log|Adj',
             caption = "Binomial (logit) regression models predicting M5S vote using survey responses") %>%
  add_header_above(c(" " = 1, "BIN (logit)" = 2)) %>%
  kable_styling(font_size = 7)

```

